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The numerical renormalization group (NRG) is rephrased as a variational method with the cost 
function given by the sum of all the energies of the effective low-energy Hamiltonian. This allows to 
systematically improve the spectrum obtained by NRG through sweeping. The ensuing algorithm 
has a lot of similarities to the density matrix renormalization group (DMRG) when targeting many 
states, and this synergy of NRG and DMRG combines the best of both worlds and extends their 
applicability. We illustrate this approach with simulations of a quantum spin chain and a single 
impurity Anderson model (SIAM) where the accuracy of the effective eigenstates is greatly enhanced 
as compared to the NRG, especially in the transition to the continuum limit. 
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The density matrix renormalization group [T] 
(DMRG), devised to improve on the Wilson's numerical 
renormalization group (NRG) [2], has become the 
method of choice to simulate one-dimensional quantum 
many-body systems at zero temperature and has found 
a large number of applications in the fields of the 
condensed matter physics, quantum chemistry and 
quantum information theory where it turned out that 
DMRG is essentially equivalent to simulating quantum 
systems in terms of matrix product states (MPS) [3]. 
In all these methods, a quantum many-body state is 
represented by associating matrices to local configura- 
tions at sites and the coefficients in the expansion over 
the configuration states are given as products of the 
corresponding matrices. The DMRG can be used to 
calculate not only the ground state but also the excited 
states through the concept of targeting where these 
matrices are chosen such that they well represent many 
states at the same time. In the NRG, on the other 
hand, the low energy states of a system are expressed 
in terms of the low energy states of a smaller system 
which is a suitable description of impurity systems with 
energy scale separation. In this context, the NRG gives 
remarkably good results and has retained the position 
as a widely used impurity solver [3] whereas the DMRG 
is rarely used to calculate the excited states [5j except 
for the spectral gap. One of the reasons is that the 
NRG is less costly as it provides D excited states at the 
cost of 0(D 3 ) as compared to DMRG with targeting 
with the cost 0(D 4 ). Presently, the DMRG is mostly 
used in its time dependent form (see e.g. [5]) which is 
also true for impurity systems [71-fTT] where also other 
density-matrix related concepts are used [I2HT4] . Still, 
DMRG's remarkable ability to calculate and optimize 
many excited states has not been used in this context. 

Even as an impurity solver, the NRG has certain limi- 
tations: the hopping terms should fall off sufficiently fast 
which means lower resolution of the spectral densities at 
higher energies e.g. Hubbard bands. Otherwise it be- 



comes inaccurate for longer chains which is seen as a vio- 
lation of the Friedel sum rule [15] . It also becomes highly 
expensive for many-band impurity problems whereas an 
additional self-consistency constraint in the context of 
dynamical mean-field theory, see e.g. |16j . calls for more 
accurate impurity solvers. What makes the NRG really 
different from the DMRG is that it does not provide any 
feedback mechanism to optimize the matrices by sweep- 
ing along the chain and a method along this line has 
already been used with quantum fields |17| . In this Let- 
ter, we identify the cost function in the NRG algorithm 
and introduce a feedback mechanism by which the states 
can be optimized in a variational way under the orig- 
inal NRG cost function. We relate the NRG and the 
DMRG by identifying the common cost function in both 
approaches show that the proposed scheme improves the 
results of both methods while retains the lower, NRG- 
like, scaling of computational costs. 

Numerical renormalization group. In the context of 
the NRG, the lowest energy eigenstates of a quantum 
many-body system on n sites, H\ip a ) = E a \ip a ), are ap- 
proximated by orthonormal states = {ip^ , • ■ ■ , } 
defining an effective low energy Hamiltonian 

^ ] = E^ n W><vi n] l (i) 

a=l 

with effective energies — (tp^\H\ip]^). The set of 
effective states is obtained recursively from S^ 1 ^, 
the effective description of a system on n — 1 sites, ex- 
tended by an extra site of a local dimension d n , resulting 
in an extended Hamiltonian e C^— 1 ^* 1 - — 1 ^ 
that, in order to prevent exponential growth, has to be 
projected to a subspace of a smaller dimension D n . In 
the NRG, this subspace is spanned by the lowest energy 
eigenvectors as Hjjj = U H H^U where H*gj € C D " xD " 
is a new effective Hamiltonian and U 6 t£.( D n-id n )xDn ^ 
an isometry matrix U ff U = 1. The set of effective states 
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FIG. 1. Tensor network representation of the (weighted) cost 
function Q in the optimization of a tensor A}^ (circles) under 
the Hamiltonian matrix product operator (squares). 

{tAIi 1 ' } is given as a MPS with an external a-index [5] 

= J2 e ^ A[1]slA[2]S2 --- A[n]Sne ^--- s ") ( 2 ) 

with the right boundary vector e Q enumerating basis 
states. The set is orthonormal due to the isometry 
constraint for matrices A^ defined as [A^']]y a y = A^J S . 
The NRG projection is equivalent to requiring that the 

[nl 

sum of the new effective energies E a is minimal which, 
according to |lj), corresponds to the trace of and 
leads to a cost function 

/(A) = tr(A H Hg™tA) where A H A = 1 (3) 

which is minimized exactly by the NRG isometry U. This 
cost function is different from other DMRG-based ap- 
proaches (e.g. [H[T2]) where the cost is determined from 
the ground state alone. 

Variational optimization. The identification of the 
NRG cost function allows us to optimize not only 
the tensor associated with the last site in the NRG- 
MPS ([2| but an arbitrary tensor in the chain. The 
cost /(A^l) = Yl a Wa \H\i^) can easily be under- 
stood as a tensor network depicted on Fig. [I] where 
the a-index enumerating states \^a) is contracted with 
the a-index in the conjugate MPS (^[r'l. The cost 
can thus be minimized by varying an arbitrary isome- 
try A® in the NRG-MPS ^ which we will illustrate 
by expressing the Hamiltonian as a matrix product op- 
erator (any ID operator can be written as a MPO) 
H = £ {s .} e ■ HM& ■ • ■ HNs«eo \s' 1 ,...,s' n ){s 1 ,...,s n \ 

where tensors H^' are represented by squares on Fig. Ill 
An arbitrary isometry A^'l is extracted out of the cost 
function ^ which now takes a form 

/(A) = ^tr(LWAR^ T A ff ) for A H A = 1 (4) 

7 

where and R^ correspond to the contractions of 
the tensor networks on the left and right side of the 
chosen site j, respectively (Fig. [lj), explicitly defined as 

[L 7 ](aV)(as) = L asls i a i and [R 7 ]p>p = Rp^pi. Fur- 
thermore, a modified weighted cost function may be con- 
sidered where the effective states are weighted accord- 
ing the their importance, e.g. by a Boltzmann factor 



w a = e~^ B ° 1 (see Fig. [lj. The cost Q is minimized 
using the conjugate gradient method with a unitary con- 
straint [18] in a variational way (i.e. the sum of energies 
can only decrease) and the computational costs scale as 
0(D 3 ), same as in the original NRG. The number of opti- 
mization steps depends on the quality of the initial state 
and the desired accuracy. 

Density matrix renormalization group. A special prop- 
erty of the NRG-MPS states is a two-fold nature of 
the external a-index which acts both as an enumerat- 
ing index as well as a virtual bond in extending the 
system for a site. For a fixed system size, the exter- 
nal index can be associated with any site in the chain 
and moved along the chain by means of the singular- 
value decomposition (SVD). Interestingly, the concept of 
a moveable external index is intrinsic to the finite size 
DMRG algorithm with targeting [T]. The basic concept 
of the finite size DMRG is as follows: split the system as 
{1, • ■ • > J-l}, {j}, + {j + 2, ...,«}, find the optimal 
tensor at site j and move to the next site j + The way 
one finds the optimal tensor at site j is by considering the 
reduced density matrix (DM) for {1, . . . , j}, obtained by 
tracing the DM of the complete system ( "universe" ) over 
the environment {j + 1, . . . , n} whereas the universe can 
be either in a pure state (ground state) or in a mixed state 
of lowest energy states as is the case with the targeting. 
Formally, the eigenstates of the universe can be written 
as|* a )=E%., ) . J+1 , r) , |^-" J - 1 )|a 3 -)|a i +i)|^ +a --"'' 1 ) 
where the index a enumerates the states and 
the optimal tensor is obtained by the SVD 

as G { i S]Sj+ir)a = J2 C A \} S3 VcV C s 3+iar - Identifying 
A lj+i]s ]+ia _ aiVla ^ ar) we recover a NRG-MPS with 
the index a associated with the site j + 1 (instead of site 
n). However, this part is ignored in the DMRG since the 
tensor at site j + 1 is obtained in the subsequent step. 
Therefore, the index a is intrinsic to but hidden in the 
DMRG and is carried back and forth along the chain. 
In the DMRG, only the central sites are optimized (the 
boundary parts are treated exactly). We can however 
sweep to the very end of the chain and end up exactly 
with the NRG-MPS ji} where the last tensor A [n] car- 
ries the a index which now allows us to enlarge the chain 
for an extra site by a NRG step. Looking back, we re- 
alize that the optimal tensor G^ SjSj+ir ) a is an isometry 
matrix minimizing the same cost function ^ as in the 
NRG context, just that it represents an arbitrary pair of 
neighboring sites instead of the last site only. Hence, the 
DMRG minimizes the same cost function as the NRG. 

We also realize that the truncation in the SVD de- 
creases the accuracy of states and it is not guaranteed 
that the optimization at the next site will recover the 
same accuracy (in fact the highest accuracy is reached 
in the center of the chain) . Therefore, the set of excited 
states obtained by the DMRG can be further optimized 
by means of the variational principles proposed earlier. 
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FIG. 2. Accuracy of the low energy spectral levels for the 
tilted quantum Ising chain (I5J) on 200 sites with h x — h z — 1, 
obtained by the DMRG (redjand optimized by the variational 
method (blue), for various bond dimensions D. 



FIG. 3. Accuracy of the eigenstates for non-interacting SIAM 
chain (6} on 70 sites (N = 68) for NRG, DMRG, and the 
variational NRG, for D £ {100, 150}. We have set A = 1.7, 
=0.1,e/ = -0.1. 



Results. The variational optimization algorithm for a 
NRG-MPS set of states is put to the test on two qual- 
itatively distinct models. First, we consider a quantum 
Ising chain in a tilted magnetic field 

n— 1 n 

H = J2 + E(>w + m?) (5) 

3=1 3=1 

which is a simple example of non-integrable quantum 
spin chains. The accuracy of the approximate eigenstates 
of |5]) is quantified by a measure (H 2 ) — (H) 2 which puts 
a lower bound to the fidelity (see supplementary mate- 
rial. We consider a chain of 200 sites, obtain the excited 
states by means of the DMRG and optimize them us- 
ing the variational approach. The DMRG provides pro- 
vides highly accurate results as seen from the Fig.[2j leav- 
ing little space for improvements. However, the results 
are not optimal and can be improved by the variational 
NRG. Higher improvement is observed for smaller bond 
dimensions D which suggests a possibility of improving 
the excited states in the regime where sufficiently large 
D are not reachable due to higher entanglement, such as 
in two-dimensional or critical (see supp. material) sys- 
tems. It is also worth noting that the computational 
costs of one DMRG sweep with targeting M states scale 
as 0(nd 3 MD 3 ), compared to either 0(MD 3 + nd 3 D 3 ) 
or 0{nd 3 D 3 ) for the variational optimization where the 
a-indcx is associated with an inner or a boundary site, 
respectively. The pre- factors to 0(MD 3 ) are similar (re- 
lated to the number of Lanczos steps) . 

As the second example we consider the Single Impu- 
rity Anderson model (SIAM) which, after a logarithmic 
discretization [3] of the conductance band, is described 
by a semi-infinite linear chain where the impurity (/ CT ) is 



coupled to a chain of fermions {cj. a ) as 

JV-i 

a 3=0 

(6) 

with n a = flf a , n = n t + n±, x = \f^, and t 3 oc A~ J / 2 
(see [1] for details). The NRG Hamiltonian Hm presents 
an effective description of the Anderson problem in the 
energy interval [A Ar+1 ,A Ar ] and only the complete se- 
quence (i?o> Hi, . . .) gives the complete spectrum. As 
mentioned in the introduction, the NRG works best for 
strong scale separation, i.e. A ^> 1 which however take 
one further away from the continuum limit A — > 1 and 
yield lower resolution at higher frequencies. Contrar- 
ily, the DMRG does not rely on energy scale separa- 
tion and thus quickly falling hopping terms but is nu- 
merically more costly, scaling as 0(D 4 ) (here M = D). 
This suggests a great potential in optimizing the existing 
NRG set of states using the variational method proposed 
in this Letter where the costs only scale as 0(D 3 ), as 
in the NRG, but, like in the DMRG, no scale separa- 
tion is required. Let us first consider the exactly solv- 
able non-interacting SIAM Hamiltonian with U — with 
a discretization parameter A = 1.7. We observe from 
Fig. [3] that for a fixed bond dimension of a NRG-MPS, 
the accuracy of eigenstates obtained by the DMRG and 
the variational optimization (vNRG) is several orders of 
magnitude better than the NRG results. Qualitatively 
the same results are obtained for the absolute errors of 
energies (see supplementary material). 

We now consider the interacting SIAM Hamilto- 
nian ^ with U = — 2ef = 0.1 on 40 sites where, un- 
like the non-interacting case, we employ the symmetries 
(total particle number of either spin). We analyze the 
improvement of the first one thousand NRG states af- 
ter variational optimization (Fig. HI) in the regime where 
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FIG. 4. Accuracy of the lowest 1000 states for the interacting 
SIAM chain (6} on 40 sites (N = 38), for NRG and the varia- 
tional NRG with symmetries. Sub-sector bond dimensions are 
taken D £ {64,100,128}, truncated to {6000,8000,10000} 
states of lowest energy. We have set A = 2, £o = 0.01, 
e f = -0.05 and U = 0.1. 



NRG already produces accurate results due to scale sepa- 
ration (A = 2). Due to the symmetries, the bond dimen- 
sion D £ {64, 100, 128} (which determines the computa- 
tional complexity) now refers to the maximal number of 
states in individual subsectors whereas the total number 
of kept states is denoted by M £ {6000, 8000, 10000}. As 
shown on Fig. |4j the accuracy of eigenstates is greatly 
enhanced by the optimization. While the cost of an op- 
timization step is comparable to the cost of a NRG step, 
a smaller bond dimension is needed to represent states 
using the vNRG as compared to the NRG. This effect 
is more pronounced for small values of A while for larger 
values of A the NRG already gives highly accurate results, 
only requiring little improvement. The results are con- 
sistent also for the fidelity of excited states and the abso- 
lute errors of energies (see supplementary material). We 
stress that the logarithmic discretization is not needed 
for the variational scheme proposed in this Letter. 

Conclusion. We have proposed a variational method 
to simulate or optimize the effective low energy descrip- 
tion of one dimensional quantum systems, introducing a 
feedback mechanism to the NRG and identifying the cost 
function. Furthermore, we have expressed the DMRG 
method with targeting as a NRG method with a move- 
able external index enumerating states and shown that 
the DMRG results can be further improved by the vari- 
ational optimization. We have tested the methods on 
the quantum Ising chain in a tilted magnetic field and 
the Single Impurity Anderson model where we observed 
significant improvement of the approximate low energy 
eigenstates. The proposed technique can be directly ap- 
plied to optimize existing NRG results [IS] and could be a 
beneficial improvement of impurity solvers in the context 
of dynamical mean- field theory, see e.g. [15] . 
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Supplementary material 

In this section we sketch the algorithm proposed in the main text of the Letter, give a proof on the validity of the 
measure used to quantify the accuracy of approximate eigenstates, and present some further results supporting the 
accuracy of the proposed method. 



Variational NRG algorithm 

1. Perform n NRG steps to obtain {ipi,i = 1, . . . , M} 
(or just choose some random states). 

2. Choose a site j := n — 1, initialize i' ' and R} n \ 

3. Calculate/retrieve tensors L^p^/ and p/ . 

4. Define G q , R 9 and X and minimize /(X) 



tMxlxM 



5. A 



[XI 



(l,s),r 



6. If converged, then stop. 

7. if (J == n), update Wj and recalculate R} n \ 
8- j : ~ j ~ 1 or i : = J + 1 an d g° to 3. 

TABLE I. Sketch of the VNRG algorithm for the variational 
optimization of a NRG-MPS set of states. 



Practical Algorithm 

In Table [I] we provide a practical algorithm for the 
variational optimization of the NRG states. Let us as- 
sume we are interested in the M lowest energy states of 
a system on sites 1 , . . . , n which we write as [5] 



ivo = E< I i A[1Isl -" AlnK i' u > 1 



(A-l) 



3 = 1 



Here, the left boundary vector is a trivial vector I = 
(1, 0,0,.. .) and fi = (0, . . . , 0, 1, 0, . . .) where the one is 
at the position /i£ {1, . . . , M}. We will only consider the 
case without employing symmetries as the generalization 
is straight forward. We will use a MPO representation 
for the Hamiltonian operator, 



block R [n] 
M energies and can be written as 



which now contains the sum over 



R 



r,q.r' = $q,0&r,r' w r- 



where w r is some chosen weight function, e.g. w r = 1 
(no weighting) or w r = c~ r (position weighting) or 
w r = c-( e t— e o) (energy weighting). In case of energy 
weighting, we use the initial approximation for the en- 
ergies and re-calculate the boundary block after each 
sweep. 

For the unitary optimization we formally reshape ten- 
sors L [3] ■ R [j] , O b ' ] , and A [j] to matrices G q , R q and X 

as [G q }(i,s),(i>,s>) = J2 P L fi,i'°pJ S > [Rgkj = R i,q,j and 
[X](; s ) r = v4p r s . The idea behind this matrix reshap- 
ing is that we can now use some black-box routine which 
minimizes the following cost function 



/(X)=^tr(G[XR 9 X T ) 



(A-2) 



under the constraint that X T • X = 1. The algorithm 
which gives the best results according to our experience 
was proposed in Ref. [T%] , 

Relation to the NRG. The MPS description of our de- 
scription ( |A-1[ ) is tightly connected to the NRG algorithm 
where one diagonalizes an effective hamiltonian if/v+i 
(Eq. 42 in jj) and obtains eigenvalues E^+i and eigen- 
vectors \w) N+1 which are related to the eigenvectors of 
the smaller system (TV sites) as (Eq. 43 in [4]) 



\w) 



N+l 



^2u(w,rs)\r;s) 



N+l- 



h= £ (i|o[ 1 K^... O N^»|i)0| s ;.)( Si |. 

We assume real arithmetics, the generalization to com- 
plex arithmetics is trivial. 

In the algorithm, we employ the fundamental concepts 
of the matrix product state formalism, such as calculating 
the blocks in a recursive way, 



- [j] 



V rb'-il A \i-i]° A \j-n*' Q ]j-iW , 



1,1' ,s,s' ,p 



and similarly for the right block R^ qr ,. The starting 

block is trivial, L' ' = Iixixi- The only step which dif- 
fers from the MPS formalism is initialization of the right 



The coefficients U(w,rs) in the above relation are iden 



tical to the tensor coefficients ^j.^" 1 ' 8 used in JA-1 1. 



Relation to DMRG with targeting. The main differ- 
ence between the DMRG with targeting (DMRGt) and 
the VNRG is that the former optimizes the set of states 
by moving the external bond along the chain whereas in 
the latter the external bond is fixed to a chosen site. In 
the DMRGt, the energies (and thus our cost function) do 
not monotonically decrease when moving the bond along 
the chain but there exists a site in the chain where the 
energies are minimal. Typically, in models with trans- 
lation invariance such as the Ising model or the Hciscn- 
berg model, the highest accuracy is reached when at the 
center of the chain. In impurity models with falling hop- 
ping terms, however, the optimal site to put the external 
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leg is (for long chains) the end of the chain - in which 
case the cheaper variant (scaling as D 3 ) of the VNRG is 
preferable. We have described the cheaper variant of the 
VNRG in the manuscript but we have skipped the more 
expensive one where the external leg is associated with 
some inner site. In fact, to reduce computational com- 
plexity, we do not actually associate the external leg to 
a site but rather introduce an additional 3-dim tensor on 
the bond connecting two sites: the third leg (the external 
one, similarly to the "physical" leg in the matrix product 
state) plays the role of the external leg. This way, the op- 
timization of the more expensive variant of the VNRG is 
as follows. All tensors associated to physical sites are op- 
timized as in the cheaper VNRG, that is minimizing the 
cost function ( |A-2[ ) under a unitary constraint. This in 
total costs nO(D 6 ) for one sweep along the chain. The 
additional tensor which selects the eigenstates is opti- 
mized in the same way as all tensors in the DMRG with 
targeting. The DMRG with targeting involves finding M 
lowest eigenstates of a sparse (Dd) 2 x (Dd) 2 matrix and 
similarly in our case where we seek M lowest eigenstates 
of a sparse D 2 x D 2 matrix. The cost of the DMRGt is 
thus nO(M(Dd) 3 ) for one sweep whereas the equivalent 
step contributes 0(MD 3 ) to the VNRG, just once m a 
sweep. 
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FIG. 5. Fidelity of the first 1000 approximate eigenstates 
of the SIAM on 40 (N = 38) sites obtained by NRG and 
vNRG and total number of kept states M € {4000,6000}. 
We have set U = -2e/ = 0.1, £ = 0.01, and A = 2 and 
used (sub-sector) bond dimension D — 5000. The basis for 
comparison are the NRG results with M = 16000 kept states 
(and D = 6000). 



The second term on r.h.s is upper bounded by e 4 



e 2 (l 



< 



G 2 



(A-7) 



On accuracy of the eigenstates 

Statement. Let H be a real symmetric linear opera- 
tor with eigenvectors {%} such that H\^j) = Ej\^j), 
and let be an approximation to ^ j with fidelity 
F = \(* j \* J )\. UF> ^,then 



1 — F < 



V2 



( th\H 2 \h) - (h\H\h 



A. V~ 3 



1/2 



(A-3) 



where Aj = naxij>^j\Ej — Ej>\ is the spectral gap to the 
nearest eigenlevel. We have assumed that all considered 
states are normalized with respect to the Euclid norm. 
We assume real algebra (generalization to complex alge- 
bra is trivial). 

Proof. Let us write the approximation VPj as 



l^yi^l^+el^) (A-4) 

where e > 0, ($|#,-) = and ($|$) = = 1. It 

is easy to show that the measure G 2 = {^fj\H 2 \^j) — 

~ 2 

(tyj\H\i£j) can be written as 

G 2 =e 2 ^ j \{H-E j f\^)-e\{^\H\^ j )-E j f (A-5) 
or as an expression for e as 

G 2 , A^H^ii-Etf 



(^\(H - Ej) 2 ^) 



(A-6) 



The denominator is lower bounded by A 2 which, assum 



ing e 2 < | and thus (1 — e 2 ) 1 < 2, finally gives us 
the upper bound on e and the lower bound on fidelity 



F = Al^e 2 > 1 - e (for e < 1) 



as 



e 2 < 2 



G 2 



F>l-V2 



G 



(A- 



The only assumption we made was e 2 < | and Aj ^ 0. 

On Figure [5] we present the fidelity for the interacting 
single impurity Anderson model (SIAM) with respect to 
very precise results obtained with the total number of 
states M = 16000 and a bond dimension D = 6000 in in- 
dividual subsectors, such that the actual number of states 
in individual subsectors after truncation to in total M 
states is always smaller than D and thus no additional 
truncation is made in subsectors. We compare the results 
obtained with M = 4000 and M = 6000, in both cases 
again using a sufficiently large D. Let us stress that we 
compare the NRG and vNRG states to the very precise 
states obtained by the NRG which were not further opti- 
mized in order to avoid possibly unfair comparison with 
preference for vNRG (no change in the fidelity is observed 
even if we do optimize the results for M = 16000). We 
observe that the accuracy of the approximate eigenstates 
is always improved by the variational method proposed in 
the manuscript. The improvement is most significant for 
those states with a low energy, since the higher excited 
states are more entangled(see e.g. V. Alba, M. Fagotti, P. 
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FIG. 6. Absolute difference of energies of the non-interacting 
SIAM (top) and the interacting SIAM (bottom) . Compare to 
Figs. 3 and 4 in the Letter. 



FIG. 7. Absolute error of energies (top) and the measure 
(H 2 ) — {H} 2 for the quantum Ising chain of 200 sites in critical 
magnetic field. 



Calabrese, J.Stat.Mech.0910, P10020, 2009 for an analy- 
sis on integrable spin systems) and require a larger bond 
dimension in the description with matrix product states. 

A single approximate eigenvalue in general does not 
tell anything about the quality of the corresponding ap- 
proximate eigenvector except for the ground state (or the 
state on the opposite side of the spectrum) . We consider 
the first M lowest lying states, so the situation might 
be similar also in our case but we are not able to prove 
it. On the other hand, good eigenvectors automatically 
also mean good eigenvalues, seen from the fact that the 
eigenvalue error is bounded from above by 

«#|iT|#) -Ef < e 4 A 2 . (A-9) 

We observe that the results for the difference in energies, 
\{H) —E\, and the measure (H 2 ) — (H) 2 are qualitatively 
very similar, as shown on Fig. [6] where we plot the error 
of energies for the same data used as in Figs. 3 and 4 in 
the Letter. 

Some additional results 

We also test the variational NRG method on an exam- 
ple of quantum Ising chain in critical transverse magnetic 



field. Similar to the tilted Ising presented in the Letter, 
we first calculate D excited states with a bond dimension 
D by means of the DMRG with targeting and then im- 
prove the results using the variational scheme. The quan- 
tum Ising chain in transverse field is a quadratic model 
and as such exactly solvable which allows us to compare 
the obtained energies (eigenvalues) to the exact values. 
We observe (Fig. [7]) a very good qualitative agreement 
between the measure (H 2 ) — (H) 2 which quantifies the 
accuracy of the eigenstates, and the absolute errors of the 
energies. Note that in this case, the improvement on the 
DMRG is more significant, since the model is critical and 
the excited states carry more entanglement - and thus de- 
mand a higher bond dimension in the MPS description. 



Expectation values 

The accuracy of the expectation values is connected 
to the fidelity. The zero temperature Green function is 
determined by matrix elements (^j ;\d a \0) and the cor- 
responding energy differences Ej — E Q where E is the 
energy of the ground state, and also |0) and the 

corresponding energy differences, see [1] for details. To 



vnrg/ 10000 


nrg/4000 vnrg/4000 


nrg/6000 vnrg/6000 


0.0122411109 


0.012240 0.0122411 


0.0122410278 0.0122411168 


0.0115338997 


0.0115338108 0.0115339045 


0.0115338277 0.0115339101 


0.0106204130 


0.0106203467 0.0106204772 


0.0106203414 0.0106204360 


0.0088452899 


0.0088452006 0.0088453015 


0.0088452206 0.0088452934 


0.0071215907 


0.0071215287 0.0071215980 


0.0071215428 0.0071215907 


0.0066356189 


0.0066202112 0.0066266177 


0.0066353419 0.0066353341 


0.0051098361 


0.0051098100 0.0051098444 


0.0051097957 0.0051098185 


0.0049692462 


0.0049637808 0.0049721638 


0.0049693508 0.0049692198 


0.0039662029 


0.0039661700 0.0039662269 


0.0039661751 0.0039662175 


0.0034311308 


0.0033144374 0.0033174038 


0.0034307370 0.0034306428 



TABLE II. The largest ten values of |(i/)j|d CT |0>| for the SIAM 
on 40 sites, U = 0.1, A = 2, £o = 0.01, e f = -0.05. The 
unreliable digits are colored red. 
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FIG. 8. Absolute error of energies in the subsector d a \0) for 
SIAM on 40 sites, U = 0.1, A = 2, £ = 0.01, e/ = -0.05. 



obtain the complete Green function, one would have 
to merge contributions from all NRG patches (i.e. for 
all system sizes) which is beyond the scope of this 
manuscript. We instead only calculate the ten largest 
matrix elements mentioned before and compare the ac- 
curacy before and after the variational optimization and 
show the results in Table [Til We observe that by the varia- 
tional optimization we gain roughly one digit of accuracy 
in the matrix elements. 

As seen from the matrix elements, the only states that 
contribute to the Green functions are those which have 
one particles less or more than the ground state, that is 
the states which live in the subsectors d a \0) and d£|0), 
respectively. We may therefore restrict the optimization 
to only these two sector which enhances the accuracy and 
reduces the computational costs. The ground state can 
be optimized independently from that by the DMRG or 
MPS techniques. In Fig. [8] we show that improvement of 
the energies in the subsector containing the state d a \0) 
where the basis of comparison are the NRG results with 
M = 16000. In this case we use sufficiently large sub- 
sector bond dimension D = 6000 such that no additional 
truncation is done in the sub-sectors and we take the 
NRG results where in total M <E {4000, 6000} states were 
kept. 

For finite-temperature Green functions we also need 
excited states in other symmetry subsectors with respect 
to the total particle number and total S z . However, they 
can all be optimized independently and in parallel. 



